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>--^ . The molasses tail in dense hard core fluids is investigated by extensive event-driven 

^N ■ molecular dynamics simulation through the orientational autocorrelation functions. Near 

^.,^' the fluid-solid phase transition, there exist three regimes in the relaxation of the pair orien- 

ts , tational autocorrelation function, namely the kinetic, molasses (stretched exponential), and 

diffusional power decay. The density dependence of both the molasses and diffusional power 
regimes are evaluated and the latter compares with theoretical predictions in three dimen- 

r^ ' sions. The largest cluster at the freezing density of only a few sphere diameter in size persist 

for only about 30 picoseconds (~ 2.8 x 10~^^[s]). The most striking observation through the 
bond orientatinal order parameter is the dramatic increase of the cluster size as the freezing 

i-^ ' density is approached. 

o 

^ . 

§1. Introduction 

2 ' The long slow decaying potential part of the shear-stress autocorrelation function 

"^ ■ (SACF) has been called the "molasses tail" to differentiate it from the hydrodynamic 

origin of the long time tail in the velocity autocorrelation function^ '^ and to empha- 
size its relation to the highly viscous glassy state.'^l''™ The decay of SACFs have been 
well investigated by mode coupling theory (MCT) and kinetic theory, both leading 

^5 \ to the same power form as in the case of the velocity autocorrelation function. This 

Q ' applies to the kinetic part of the shear autocorrelation function, but not to the po- 

CJ . tential part which is central to the molasses tail. From the numerical point of view, 

the long time tail of SACF has been computed by several molecular dynamics (MD) 
simulations in the early 80s; Evans (1980),^' Wood & Erpenbeck(1981),^'and Morriss 
& Evans(1985).'^ These studies seem to show that the long time tail of both the 

^3 ! kinetic and potential parts have a power decay consistent with MCT, however, the 

''si" ' amplitude of SACF in dense fluids was found to be orders of magnitude greater than 

predicted. This discrepancy was ascribed to the possibility that the numerical results 

Ji^ \ were not run long enough. However, we show that the tail is caused by slow struc- 

(<— ^ ' tural relaxation in the dense liquid around the peak of the structure factor rather 

than by hydrodynamic phenomena at long wave length. This was the starting point 
of the theory of glassy transition.'^''^ Since then many papers focusing on density 
fluctuations in the study of the glass transition have been published, however, the 
microscopic mechanism of the stress field relaxation have not been clarified yet. 

5^ \ Twenty years ago, Ladd and Alder have speculated that the long time tail of the 

shear stress auto-correlation function near the solid-fluid transition point in the hard 
sphere system is due to transient crystal nuclei formation. ^'^' They found that the 
potential part of the SACF and the angular orientational auto-correlation function 
(OACF) are identical in the long time limit and show non-algebraic decay in time. 
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Since the evidence suggested that the reason for non-algebraic decay is structural 
relaxation rather than hydrodynamic flow, an attempt was made to understand this 
slow decaying mechanism by decomposing the OACFs into two-, three-, and four- 
body correlations, however, many of these correlation functions, especially the four- 
body correlations, have not been obtained accurately due to the computer limitation 
at that time. Even with today's much better computers, it is not possible to get 
accurate information on the four body distribution and hence we made an attempt 
to use what could be learned about the cluster size from the bond orientational 
parameter. The results, presented here, look promising, and we plan to extend, this 
nearest neighbour bond order parameter to further neighbours. The prediction of 
the cooling rate necessary to prevent crystallization requires knowledge of the rate 
of growth of a cluster the size of a critical solid nucleus and the time it exists in 
this transient state. Therefore, analysis of this slow decaying process of the OACF 
and its decomposition into several components is expected to be a key factor in 
understanding the onset of the glass transition.'^''™ 

Previously,^'' we have reported on a two dimensional system consisting of elastic 
hard disks at a single density near the solid-fluid transition point placed in a square 
box with periodic boundary conditions, using a modern fast algorithm based on 
event-driven MD simulation.li^I We confirmed the non-algebraic decay (stretched 
exponential) at intermediate times presumably due to the existence of various sized 
solid clusters at high densities decaying at different rates. We also determined the 
length of time for which the biggest such nuclei exists. In this study, we focus on 
the rapidly increasing time with increasing density for the decay of OACFs and are 
able to establish the length of time for which the biggest such nuclei exists at each 
density. We also compare the results to theoretical predictions in the subsequent 
power law decay. We will further report results for hard spheres for different particle 
numbers. We use a more efficient program code for calculating pair contributions to 
the OACFs, allowing for more accurate results. 

§2. Decomposition of the Orientational Autocorrelation Function 

We focus on the potential part of the SACF (J^(t) J^(0)), where J^ is the 
potential part of the momentum current J^, where the molasses tail appears. Since 
velocities and positions are no longer correlated beyond a few mean collision times, 
only the orientational part of SACF, namely the OACF {Oxy{t)Oxy{0)) needs to be 
studied.i°''ii' O^yit) is defined as 

0.y{t) = Y,^6{t-Q, (2-1) 

c 

where Ylc n^sans contribution at the colliding time tc at which {xij,yij) = {xi — 
XjiVi — Vj) are the relative positions between hard spheres (disks) i and j. 

To avoid the delta function singularity of Oxyif) for hard particles, the alterna- 
tive Einstein-Helfand expression^ 'l^ involving the second derivative, obtained by 
the numerical differentiation with five point stencils, is adopted for calculating the 
correlation function. 
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2dt^ 



Cit) = {0,y{t)0,ym = --^{(Git) - cm'), (2-2) 



where 



Git) = Y,^0{t-t,), (2-3) 

c 

and where 0{t) is the unit step function. There are three independent orientational 
factors {Oxy, Oyz, and Ozx) in 3D. We introduce the reduced (e.g. Enskog values) 
orientational function C*{s) in terms of the reduced time s = t/tQ, 

C*is) = ^C(t), (24) 

Vr]E 

where to is the mean free time and rjE is the Enskog shear viscosity. For a hard 
core fluid, the Enskog shear viscosity can numerically be estimated by the expres- 
sion of Gass(1971j^^' with the third Sonine polynomial approximation in 2D and 
Wainwright(1964)^^' with the equation of state by Ree&:Hoover(1967/^ in 3D. 

The total correlation function C{t) can be decomposed into pair C2{t){ij — ij 
collision), triplet C3(t) {ij — ik collision), and quadruplet Ci{t) {ij — kl collision) 
contributions,'^ 'tl^ where i,j,k,l are particle index. 

C{t) = C2{t) + C3{t) + C4t). (2-5) 

The pair contribution C2{t) is defined as, 

INN \ , r2 / N N \ 

C2it) = J^{i: E 0%{t)0%{0)\ = —^/Y: E GHt?), (2-6) 

\ i j(j>i) I \ i j{j>i) I 

since G*-'(0) = {G^^{t) = Oxy{t)). To calculate C2, we introduce a "collision pair 
index" c^ = (q — 1)A^ — Cj(cj — l)/2+Cj— q, where A^ is the number of particles, which 
identifies a given pair quickly, avoiding having to check whether the collision pair is 
same as before in the process of calculating G^^{t). This speeds up the calculation 
considerably. O 

§3. Results 

The system consists of hard disks or hard spheres placed in a L^ x L^, (= A) square 
box or a Lx X Ly X Lz{= V) cubic box, respectively, with periodic boundary condi- 
tions. Initially, the simulation systems for each packing fraction v {= N7r{a/2)'^ /A 



*' For example, in case of A'^ = 4, the total number of collision pairs are N{N — l)/2 — 6, 
which can be listed as {ci,Cj) — (1,2), (1,3), (1,4), (2,3), (2,4), (3,4), where d < Cj. By using the 
"collision pair index", we obtain Ck = 1,2,3,4,5,6 for each collision pair, respectively. Therefore, 
it is easy and convenient to deal with the collision pair as the sequential number to sort and insert 
the array of correlation pair G'-' (i) when we make the computer program. 
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for 2D, or (4/3)A^7r(o"/2)^/y for 3D), are prepared as the equilibrium state by a 
sufficiently long preliminary run. Relative to the close packed area Aq, and vol- 
ume Vo, V = vr/(2\/3(A/74o)) and u = 7r/(3\/2(V^/Vo)), respectively. The system 
evolves through collisions, using an algorithm based on event-driven Molecular Dy- 
namics (MD) simulation. ^^' Most of the calculations are done with the particle 
number N = 512 and 4, 096. The density is varied at relatively dense value near the 
solid-fluid transition point, which are known to be v^ ~ 0.70 for 2D and ~ 0.48 
for 3D. Ladd & Alder (1989)1^ have calculated OACFs and its separate parts 
at several densities; A/Aq = 1.35,1.4,1.6 for 2D and V/Vq = 1.5,1.6,1.8,2.0 for 
3D. These densities correspond to i^ = 0.672 • • • , 0.648 • • • , 0.567 • • • for 2D and 
V = 0.493 • • • , 0.463 • • • , 0.411 • • • , 0.37 • • • for 3D, respectively. In our calcula- 
tion, we use several additional packing fractions up to solid-fluid transition density; 
u = 0.69, 0.67, 0.65, 0.57 for 2D and v = 0.47, 0.45, 0.40, 0.35 for 3D. 

3.1. Density & system size dependence of OACFs in 2D and 3D 

First, we calculate the total OACFs for different packing fractions and system 
size to investigate the functional form of the molasses tail at long times, and the 
effects of periodic boundary conditions. In Figs. 1 and 2, the dimensionless OACFs 
for various packing fractions near the solid-fluid transition point in terms of reduced 
time s = t/tQ for 2D and 3D are shown, respectively. Statistical averages are made 
over 10^^ to 10^^ total collisions of event-driven MD runs, which are 10^ ~ 10^ larger 
than previous work.l^ Typical error bars are also shown at long times. Two different 
number of particles are shown in both 2D and 3D to observe periodic boundary effect 
due to sound wave propagation across the system. The sound wave propagation time 
explicitly appears in the long time tail of velocity auto-correlation functions (e.g., 
Ref.^), but appears to have no effects here, as expected, since OACF does not involve 
velocities. 

The OACFs for both 2D and 3D clearly show the slower decay when the system 
become denser and their decay are not the power law form as for the hydrodynamic 
tail ~ s^ ' ^, where d is dimension. These data suggest that the theoretical prediction 
of MCT for the long time tail must be reconsidered in dense liquids. Based on these 
results, we choose the relatively smaller system size N = 4096 in 2D and N = 512 
in 3D for efficiency to calculate C2 in the following subsection. 

3.2. Density dependence of C2 in 2D and 3D 

To investigate the density dependence of C2 up to the time of the diffusional 
power regime, long runs were performed for various packing fractions. Since the 
maximum collision index (i.e., Cfc = N{N — l)/2) increases as ~ 0{N'^), it is difficult 
to calculate correlation functions averaged for all possible collision pairs c^ due to 
the restriction on CPU time and memory, if N becomes more than one thousand 
particles. For instance, one needs about 100M6ytes of memory for the calculation of 
C2 for all possible collision pairs in the 500 particles system, comparable to the Ladd 
& Alder(1989) calculation, which presented no serious problem. However, in the case 
of A^ = 4, 096 the maximum collision index becomes c^ = 8, 386, 560, which needs 
about WGbytes of main memory. To reduce the calculation to a manageable one, we 
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Fig. 1. A log-log plot of the dimensionless OACFs (see Eq.(2.4)) for various packing fractions for 
two different sized systems in 2D versus reduced time, s = t/to are shown. Arrows describe the 
approximate values of the sound wave transversal time (the left for N = 1024 and the right for 
N = 4096.) 
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Fig. 2. As Fig.l except for 3D. 
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restrict the range of collision pairs Ck from 1 to lOA^ in the averaging procedure. This 
restriction of c^ only affects to error bars, but we perform sufficiently long time MD 
runs for C2 , much longer than the typical time scale of the largest cluster breaking 
up. 




Fig. 3. CJs for various packing fractions near the solid-fluid transition point in 2D (left) and 3D 
(right) in terms of reduced time s = i/io- Arrows indicate the sound wave transversal time 
for each system size. Theoretical slopes by using the effective diffusion constant sd v = 0.57 in 
2D(left) and the diffusion constant obtained from actual simulation for each packing fraction in 
3D(right) at long times are given by the lines (see, Sec. 3.4). 

We performed up to 10^^ collision for several independent runs for different time 
resolution, namely At = to, lOto, and lOOto at each packing fraction. We confirmed 
that these runs give consistent results. The efficiency of these calculation in such that 
it does not significantly slow down the MD calculation. The error bars in the figures 
are smaller than the symbols except at a few points at long times. As indicated in 
Fig. 3, C2 in both 2D and 3D decays slower as the density becomes higher. The 
results are consistent with both our previous method at v = 0.69^ and Ladd & 
Alder (1989).'^ It is interesting to observe that qualitatively the decay is similar for 
2D and 3D. Next, we investigate the decay of C2 in the two separate time regime, 
namely the molasses and diffusional regimes. 

3.3. Molasses regime in C2 

We investigate the scaling behavior in the molasses regime as we did previ- 
ously.^ In the molasses regime, C2 was found to have a scaling law of the stretched 
exponential form as. 



Cl{t) = 5exp 



(3-1) 



where ^2 is the relaxing time in the molasses regime and a is a density "independent" 
stretched exponent. The value of the exponent a was found to be ~ 0.80 for 2D 
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Fig. 4. Decay of the CI in a semilog plot for 4096 particles against the scaled time {t/t2f'^ for the 
several packing fractions in 2D. 



and ~ 0.70 for 3D in previous workflSli By adopting a simple fitting method to C2 
as given by eq. (|3-1|) with the given values of a, we estimated ^2- C2 for the several 
packing fractions plotted against the scaled time (t/t2)" are shown in Figs. 4 and 5. 
We found that the curves for various packing fraction collapse to a straight line in 
the regime of (i/i2)" between 1 and 5. Thus, the decay of C2 is in good agreement 
with the stretched exponential form of eq. ()3-ip . which is a similar result given in 
Ref.l^ The crossover time tcs from the molasses regime to the diffusional power 
regime in terms of the mean free time io are summarized in Table I and Fig. 6. 



packing fraction(2D) u 0.57 0.65 0.67 0.69 

crossover time tcs /to 100 287 536 1364 

packing fraction(3D) u 0^35 0^40 045 0.47 

crossover time tcs/to 55 115 281 431 

Table I. The crossover time tcs from stretched exponential to power form in terms of the mean 
free time to for various packing fractions are shown. 

The log of the crossover time is increasing linearly in 3D, but faster than exponential 
in 2D with packing fraction near the solid-fluid transition point Vc- 

3.4. Diffusional regime in C2 

In the regime of (t/i2)° > 6, the data show a different functional form from 
the molasses regime, in which the correlation function changes from a stretched 
exponential to a power decay. 



cm = B't-^, 



(3-2) 
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Fig. 5. Decay of the C2 in a semilog plot for 512 particles against the scaled time (t/t2)'^'^ for the 
several packing fractions in 3D. 
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Fig. 6. The log of the crossover time tea divided by the mean free time to and the log of the 
amplitude B' , both as a function of the packing fraction given relative to the freezing point. 



where /3 and B' are the exponent and the amplitude, respectively. From past the- 
oretical work based on MCT and the pair diffusion mechanism, it is predicted that 
/3 = 7/2 in 3D. Leegwater&:Beijeren(1989p^ investigated the equation for two par- 
ticles diffusing relative to each other, which predicts at long times the behavior of 
C2 as 
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(3-3) 

(34) 



where n = N/V is the number density, D is diffusion constant, E stand for En- 
skog, g{a) is the pair distribution function at contact, and 1^/2 is a modified Bessel 
function. They also found that a modified theoretical prediction by incorporating 
recollision contribution and a potential of mean force is in excellent agreement with 
the MD data throughout the intermediate-time and long-time regimes t/to > 10, 
although this requires an effective diffusion constant about 10% less than for a single 
particle.'^ By the same theoretical argument applied to the 2D, '221' -we obtain 



6APr,E \4DJ ^ '' 

In 2D, the theoretical decay is C2{t) ~ t~^ in the long time limit. 
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Fig. 7. Decay of CJ in the diffusional power decay regime in terms of (i/i2)~ for the several 
packing fractions in 2D(left) and 3D(right) for 4096 and 512 particles are shown. The solid lines 
are theoretical results for v — 0.57 in 2D(left) and v — 0.47 in 3D(right), respectively. The dot 
lines are the linear fitting lines. 



Plotting C| in terms of (t/t2 



-7/2 



in Fig. 7, in the 3D case the data lead to a linear 



function, which indicates consistency only with the exponent of theoretical prediction 
of /3 = 7/2 in the long time limit. In the 2D case, numerical data suggest the exponent 
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deviates from /3 = 3. Because of the divergence of the diffusion coefficient in 2D, 
the entire theory needs to be reworked. The amphtudes of the simulation data in 
3D are fairly good agreement with the theoretical results incorporating recollision 
contributions and a potential of mean force through the extra factor of 9(0'), (i.e., 
(25/9)g((T)C|(i) in the long time limit), as shown in Fig. 3, 6 and 7, which are 
described in Ref.^^' In Table II, the values of amplitude B' of eq. (3.2) are shown. 
We also show B' in terms of scaled packing fraction {v/vc- ^c is solid-fluid transition 
point) in Fig. 6. 



packing fraction(2D) v 0.57 0.65 0.67 0.69 

Amplitude B' 1.42 17.8 75.5 1849 

packing fraction(3D) v 035 040 0^45 0.47 

Amplitude B' 0.12 0.47 3.16 8.57 

Amplitude(Tlieory) B' 0.047 0.32 4.00 11.58 

Table II. The amplitude B' for various packing fractions are given and plotted in Fig. 6. 



§4. Discussion 

The cause of the stretched exponential relaxation in the molasses regime is con- 
sidered to be due to the distribution of different life times of transient clusters of 
nuclei in dense fluid systems. Thus, there exist several exponential relaxation decays 
for each collision pair. To confirm the above speculation for the origin of the molasses 
effect in 2D at the microscopic level, we calculated the bond orientational order, ^% 
as an alternative to C2, and furthermore, to visualize the distribution of crystal clus- 
ters in the hard disk system. The usual bond orientational order parameter ^^ for 
each hard disk i is defined by 



1 ^' 






(4-1) 



where Ni is the number of the nearest neighbors around the tagged particles i, and 
6ij is the angle between the position vector from the disk j to i and an arbitrary fixed 
reference axis (e.g., x-axis). Two disks are defined as neighbors if the distance is less 
than y/Sa. In this definition, (J)q takes on values between and 1, which measures 
the degree of crystallization in terms of considering only nearest neighbors. 

Figure 8 shows typical snapshots of the spatial distribution of (pQ. The gradation 
in shading of the particles indicate the value of (J)q, the darker, the closes to unity. 
We clearly observe the dramatic growth of several solid nuclei as the density nears 
solidification. The averaged values of (p^ and the number of particle with 0.9 < 4>q < 1 
for each packing fraction are summarized in Table III. 

Understanding the slow decaying process of the OACF is considered a key factor 
in understanding the onset of the glass transition, and is here qualitatively confirmed 
to be due to the transitory existence of solid nuclei in 2D in Fig. 8. Recently, their 
existence was also demonstrated in colloid experiments .1^ For quantitative purposes 
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Fig. 8. The spatial distribution of (Jjq for 4096 particles system at a given time for two packing 
fractions, u — 0.65(left) and 0.69(right). The darker the region, the closer (pa is to unity. 



packing fraction(2D) 


V 


0.57 


0.65 


0.67 


0.69 


averaged bond orientational order 


4>l 


0.36 


0.42 


0.44 


0.47 


particle number with 0.9 < ^g < 1 


Ns 


11 


58 


82 


131 



Table III. The averaged bond orientational order <j!>g and the number of particle with 0.9 < (f)\, < 1 
for various packing fractions are shown. 



it is crucial to determine the spatial extents of the crystal nuclei and how long they 
persist. In that direction we quantitatively determined in this paper that the decay 
of the pair distribution of C2 could be characterized by a relaxation time ^2 at 
intermediate time (molasses regime) and later on by a pair diffusional breaking up 
of the clusters. At the quantitative crossover time shown explicitly in Table I and 
Fig. 6, the largest cluster at the freezing density of only a few sphere diameter in 
size persist for ~ 431to) which corresponds for typical argon parameters estimate of 
to ~ 6.4 X 10~^^ [s] to only about 30 picoseconds (~ 2.8 x 10~^^[s]) in real liquids. 
To make the further quantitative progress, we need to investigate the orientation 
correlation function of the quadruplet component, C4, as a function of the distance 
between the two colliding pairs. When this is carried out as a function of density, it 
could be possible to tell how the cluster size distribution changes and how long they 
last as the solidification density is approached, and, subsequently, determine how fast 
one has to increase the density to get a glass instead of a crystal. This is however, 
computationally a very demanding task. Instead we are planing to investigate an 
extension of (/)g to a higher order orientational parameter involving further neighbors 
and thus, hopefully, improve the methodology of an alternative efficient calculation 
for Ci{t). 
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